source("global.R")
library(ggrepel)
library(ComplexHeatmap)
lipid <- readRDS("../data/untargeted_singletp_fin4_all.rds")HDL lipidome quantification (from untargeted)
Analysis based on the targeted lipidomic results. File: mx 733759_Agus_lipidomics_human plasma lipoproteins_07-2023 submit.xlsx. The peak heights were converted to ng/mL based on the amount and peak height of iSTD.
Data processing
Peak heights across samples were normalized using all iSTDs.
Features has been filtered by WCMC regarding the peak heights in samples and blanks.
Features with InChI Keys that couldn’t be found in the PubChem were excluded.
Compounds with identical InChI Keys were selected by lowest CV in pooled samples.
Assign lipid class from LIPID MAPS and manually curated table.
Concentration were calculated based the istd with most similar unsaturation and then chain length.
Cholesterol istd is missing. Import the value from mx733759_Agus_Lipids_Single-point quant_with cholesterol_Submit_10-17-23.xlsx
Total peak heights consists of major lipid species which can be quantified with iSTDs and minor species which don’t have corresponding iSTDs. The percentage of both types of species are shown below.
peak <- lipid$peak
lipid <- lipid$quant
peak$fdata$major <- ifelse(featureNames(peak) %in% featureNames(lipid), "Major", "Minor")
lipid.major <- summarize_feature(peak, "major")$edata
lipid.major <- apply(lipid.major, 2, function(x)x/sum(x)*100)
df <- as.data.frame(lipid.major) %>%
rownames_to_column("Type") %>%
pivot_longer(cols = contains("_0"), names_to = "fraction", values_to = "Perc") %>%
mutate(fraction = factor(fraction, levels = colnames(lipid.major))) %>%
group_by(fraction) %>%
mutate(ypos = 100 - cumsum(Perc) + 0.5 * Perc)
ggplot(df, aes(x = "", y = Perc, fill = Type)) +
geom_bar(stat = "identity", color = "white") +
geom_label(aes(y = ypos, label = round(Perc, 1)), fill = "white") +
coord_polar("y", start = 0) +
theme_void() +
scale_fill_npg() +
facet_wrap(~fraction, nrow = 3)
The quantified data only included major lipid species.
CVs of pooled samples
The CVs of pooled samples for each lipid species were calculated using the following formula:
\[ CV_{lipid\ species} = SD_{lipid\ species}/Mean_{lipid\ species} * 100\% \]
We are less confident with the lipid species with a high CV. These lipid species will be excluded for further analysis.
hist(lipid$fdata$qc_cv, main = "The distribution of CVs")
sum(lipid$fdata$qc_cv>30)[1] 11We excluded lipid species with a CV > 30%. 11 features were excluded.
lipid <- subset_features(lipid, lipid$fdata$qc_cv<30)The hierarchical relationship of lipid classes.
Following analysis included all quantified lipid species.
Included all quantified lipid species. 456 lipid species for analysis.
Lipid class
Pie chart
lipid.class3 <- subset_features(lipid, !is.na(lipid$fdata$
class))
lipid.class3 <- summarize_feature(lipid.class3, "class")
edata.prop3 <- apply(lipid.class3$edata, 2, function(col){
col/sum(col)*100
})
# pie chart
pies.class3 <- lapply(sampleNames(lipid), plotPie, edata.prop3)
legend.class3 <- get_legend(pies.class3[[1]]+theme(legend.position = "bottom"))
pies.class3 <- lapply(pies.class3, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class3, nrow = 3)
Bar plot
# bar plot
plotBar(edata.prop3)
PCA
edata.scaled3 <- scale(t(edata.prop3), center = T)
res.pca <- PCA(edata.scaled3, scale.unit = T, graph = FALSE)
# fviz_eig(res.pca, addlabels = TRUE)
fviz_pca_ind(res.pca,
axes = c(1,2),
repel = TRUE,
title = NULL
) +
theme_cynthia_bw() +
theme(
legend.position = "none"
) +
labs(title = NULL)
Lipid main class
Pie chart
lipid.class2 <- subset_features(lipid, !is.na(lipid$fdata$Main.class))
lipid.class2 <- summarize_feature(lipid.class2, "Main.class")
edata.prop2 <- apply(lipid.class2$edata, 2, function(col){
col/sum(col)*100
})
# pie chart
pies.class2 <- lapply(sampleNames(lipid), plotPie, edata.prop2)
legend.class2 <- get_legend(pies.class2[[1]])
pies.class2 <- lapply(pies.class2, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class2, nrow = 3)
Bar plot
# bar plot
plotBar(edata.prop2)
Lipid category
Pie chart
# edata
lipid.class <- subset_features(lipid, !is.na(lipid$fdata$Categories))
lipid.class <- summarize_feature(lipid.class, "Categories")
edata.prop <- apply(lipid.class$edata, 2, function(col){
col/sum(col)*100
})
# pie chart
pies.class1 <- lapply(sampleNames(lipid), plotPie, edata.prop)
legend.class1 <- get_legend(pies.class1[[1]])
pies.class1 <- lapply(pies.class1, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class1, nrow = 3) 
Bar plot
# bar plot
plotBar(edata.prop)
Following analysis only included lipids curated in LIPID MAPS.
220 lipid species for the analysis.
LM: sub class
Pie chart
lipid.class.sub <- subset_features(lipid, !is.na(lipid$fdata$LM_sub_class))
lipid.class.sub <- summarize_feature(lipid.class.sub, "LM_sub_class")
edata.prop.sub <- apply(lipid.class.sub$edata, 2, function(col){
col/sum(col)*100
})
# pie chart
pies.class.sub <- lapply(sampleNames(lipid), plotPie, edata.prop.sub)
legend.class.sub <- get_legend(pies.class.sub[[1]]+theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 5)))
pies.class.sub <- lapply(pies.class.sub, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class.sub, nrow = 3) %>%
plot_grid(legend.class.sub, nrow = 2, rel_heights = c(3,1))
Bar plot
# bar plot
plotBar(edata.prop.sub)
PCA plot
edata.scaled.sub <- scale(t(edata.prop.sub), center = T)
res.pca <- PCA(edata.scaled.sub, scale.unit = T, graph = FALSE)
# fviz_eig(res.pca, addlabels = TRUE)
fviz_pca_ind(res.pca,
axes = c(1,2),
repel = TRUE,
title = NULL
) +
theme_cynthia_bw() +
theme(
legend.position = "none"
) +
labs(title = NULL)
LM: main class
Pie chart
lipid.class.main <- subset_features(lipid, !is.na(lipid$fdata$
LM_main_class))
lipid.class.main <- summarize_feature(lipid.class.main, "LM_main_class")
edata.prop.main <- apply(lipid.class.main$edata, 2, function(col){
col/sum(col)*100
})
# pie chart
pies.class.main <- lapply(sampleNames(lipid), plotPie, edata.prop.main)
legend.class.main <- get_legend(pies.class.main[[1]]+theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3)))
pies.class.main <- lapply(pies.class.main, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class.main, nrow = 3) %>%
plot_grid(legend.class.main, nrow = 2, rel_heights = c(4,1))
Bar plot
# bar plot
plotBar(edata.prop.main)
LM: core
Pie chart
lipid.class.core <- subset_features(lipid, !is.na(lipid$fdata$LM_core))
lipid.class.core <- summarize_feature(lipid.class.core, "LM_core")
edata.prop.core <- apply(lipid.class.core$edata, 2, function(col){
col/sum(col)*100
})
# pie chart
pies.class.core <- lapply(sampleNames(lipid), plotPie, edata.prop.core)
legend.class.core <- get_legend(pies.class.core[[1]])
pies.class.core <- lapply(pies.class.core, function(x)x+theme(legend.position = "none"))
plot_grid(plotlist = pies.class.core, nrow = 3)
Bar plot
# bar plot
plotBar(edata.prop.core)
Lipid species enriched in fractions
The analysis was based on peak heights. We assumed that the sum of peak heights in each sample represented the total lipid amount. Each lipid species was normalized to the total lipid amount before comparing across samples.
CVs of pooled samples
hist(peak$fdata$qc_cv, main = "The distribution of CVs")
sum(peak$fdata$qc_cv>30)[1] 18We excluded lipid species with a CV > 30%. 18 features were excluded.
peak <- subset_features(peak, peak$fdata$qc_cv<30)
major.classes <- c("Cer", "CAR", "CE", "Cholesterol", "DG", "FA", "LPC", "LPE", "PC", "PE", "SM", "TG")
classes <- peak$fdata$class
classes[!classes %in% major.classes] <- "Other"
peak$fdata$class <- classes
peak.rel <- apply(peak$edata, 2, function(x){
x/sum(x)
})Species level
Class level
peak.class <- summarize_feature(peak, "class")
peak.class.rel <- apply(peak.class$edata, 2, function(x){
x/sum(x)
})
p.class <- enrich_all_fraction(bg.fraction, peak.class.rel, peak.class, TRUE)
# Heatmap of enrich species in each fraction
mat.class <- as.matrix(p.class)
mat.class[mat.class<0.01] <- 0.01
mat.class <- -log(mat.class)
cols <- RColorBrewer::brewer.pal(name = "Paired", n = 12)
cols <- colorRampPalette(cols)(13)
cols <- setNames(cols, rownames(p.class))
Heatmap(mat.class, name = "-log(p)",
right_annotation = rowAnnotation(class = rownames(p.class), col = list(class = cols)),
show_row_names = FALSE)